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Abstract. A technique for the analysis of low-low intersatellite range-rate data in a gravity mapping 
mission is explored. The technique is based on standard tracking data analysis for orbit determination 
but uses a spherical coordinate representation of the 12 epoch state parameters describing the baseline 
between the two satellites. This representation of the state parameters is exploited to allow the 
intersatellite range-rate analysis to benefit from information provided by other tracking data types 
without large simultaneous multiple data type solutions. The technique appears especially valuable 
for estimating gravity from short arcs (e.g. less than 15 minutes) of data. Gravity recovery simulations 
which use short arcs are compared with those using arcs a day in length. For a high-inclination orbit, 
the short-arc analysis recovers low-order gravity coefficients remarkably well, although higher order 
terms, especially sectorial terms, are less accurate. Simulations suggest that either long or short arcs 
of GRACE data are likely to improve parts of the geopotential spectrum by orders of magnitude. 

Key words: Geopotential determination — Satellite geodesy Satellite-to-satellite tracking 
GRACE. 


1. Introduction 

For more than three decades the geodetic community has realized that satellite-to-satellite tracking 
(hereinafter SST) provides extremely strong observational constraints for determining the geopotential 
(e.g. Wolff, 1969; Vonbun, 1972). High-low satellite configurations have proven, valuable in the past 
(e.g., Kahn et ah, 1982), and they continue to do so today (Lemoine et al., 1998b; Schwintzer et 
al., 2000). Low-low satellite configurations are expected to yield orders of magnitude improvements 
in geopotential definition (National Research Council, 1997), and such a system may finally come to 
fruition in the near future with the Gravity Recovery and Climate Experiment (GRACE) mission 
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(Tapley and Reigber, 2000). 

Methods for deducing the geopotential from low-low SST data have been developed by many groups 
over the past two decades (e.g., Douglas et al., 1980; Kaula, 1983; Wagner, 1983, 1987, Colombo, 1984, 
Jekeli and Upadhyay, 1990). Many of these methods were based on semi-analytic theories that made 
various simplifying assumptions to overcome computational limitations. For example, both Kaula s 
and Colombo’s methods assumed perfectly polar orbits a “tail-biting orbit was Colombo s colorful 
phrase— which allow, among other advantages, fast Fourier techniques. Computational limitations 
are still an important consideration today, but not nearly so much as when these earlier works were 
written. 

This paper explores a direct approach to estimating gravity from SST data, relying more heavily 
on numerical integration methods than on analytic or semi-analytic methods. We present extensive 
simulations with the following SST scenario: near-polar low-low satellites, separated by a couple 
hundred kilometers, each satellite tracked continuously by GPS and each deployed with accelerometers 
to correct for non-conservative forces. The fundamental measurement is the satellite-to-satellite range 
rate, measured with a precision of order 10 -6 ms 1 ■ 

In addition to establishing a viable technique for the handling SST data, this paper addresses two 
key questions for any practical data analysis: 

1. To what extent is the SST gravity inversion insensitive to ephemeris accuracy? Specifically, is it 
sufficiently insensitive that the orbit determination and the gravity inversion can be performed 
in separate, independent steps? 

2. To what extent is the SST gravity inversion dependent on arc length? In particular, if the ac- 
celerometers are incapable of removing all non-conservative forces, including thrusting events, 
the satellite data might necessarily be broken into very short arcs. How will the gravity estima- 
tion be affected? 
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Concerning (1), our simulations described below show that, indeed, the SST data can be handled 
(within limitations) independently of the GPS data. We validate a two-step method. The first step 
concentrates only on achieving orbit accuracy (no gravity estimation). The second step uses only 
the SST data to refine certain components of the orbit while estimating gravity coefficients. A great 
advantage of this is that the orbit determination task in an SST mission can concentrate on ephemeris 
accuracy, using techniques such as empirical accelerations that are normally prohibited in standard 
satellite gravity estimation. A second practical advantage is, of course, that the gravity inversion task 
can be performed without the considerable simultaneous data processing chores associated with GPS 
orbit determination. 

The next section compares the role of arc length in conventional orbit determination based gravity 
estimation with the role it may play in a GRACE-like SST mission. Section 3 describes a transforma- 
tion that we apply to the standard orbit parameters to enable our decoupled analysis of the SST data. 
Section 4 gives the rationale behind our data simulation procedures as well as a detailed description 
of those procedures. Sections 5 and 6 present results of parameter estimations using simulated data. 
Section 5 concentrates on the estimation strategy for orbit parameters, while Section 6 presents the 
results of gravity field recoveries with orbit parameters adjusting simultaneously. 

2. Arc Length in Gravity Estimation 

When satellite tracking data are analyzed for geophysical parameter recovery, the analysis is often 
part of a simultaneous solution for orbit parameters. In such settings it is usually necessary to group 
the tracking data into “arcs” of data which span multiple revolutions of the satellite. This is because 
most tracking data types do not provide enough geometric strength to provide a unique solution for the 
orbit until a substantial portion of the trajectory is examined. Furthermore, force model parameters 
like gravity coefficients build up sensitivity with arc length, because gravity signal is usually detected 
through its effect on the trajectory of a satellite. If two trajectories for a satellite are computed, 
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each starting from the same initial conditions but using different gravity models, it usually takes 
some period of time from the initial epoch before differences in the trajectories are large enough that 
differences in the gravity signal can be inferred. The trajectories are used to compute “theoretical” 
values of the tracking data observations to which the actual observations are compared. The difference 
in the trajectories due to gravity needs to affect the computation of the theoretical tracking data 
values above the level of the precision of the actual tracking data. Furthermore, the differences in 
trajectories computed with different gravity models are often diminished because the initial state 
estimation process at some level accommodates gravity errors. Gravity is really an “indirect effect” in 
conventional tracking data analysis. Arc length is required for the indirect effect to make its presence 
felt. 

In general, long arcs are desirable. On the other hand, as the arc length grows, so does the effect of 
unmodeled forces. Therefore, the length of an arc in an orbit solution cannot be extended indefinitely 
without degrading the solution. Choice of arc length is a key decision in the analysis of tracking data. 
It depends on the geometric strength of the tracking data, the magnitude of unmodeled forces, and the 
sensitivity of the data to the geophysical parameters of interest. Gravity models like EGM96 (Lemoine 
et al., 1998a) have typically used tracking data analyzed in arcs of data from 1 to 30 days. With very 
few exceptions the tracking data used in EGM96 would not support the extraction of gravity signal 
from an arc of data significantly shorter than a day. 

In a gravity mapping mission where there is very precise tracking of the change in range between 
two satellites in the same orbit plane, many of the conditions that have always been a factor in deciding 
arc length for conventional tracking data will be much different. The most obvious difference is the high 
precision of the data. Small trajectory changes are detectable with more precise data. Furthermore, 
the intersatellite observation is, in some sense, a direct measurement of the difference in forces which 
the satellites experience at each instant, and hence almost a direct measurement of gravity (Colombo, 
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1984). The extraction of gravity signal from such a measurement is unlikely to require the type of 
trajectory analysis that conventional tracking data requires. A precise range change measurement 
between two satellites in the same orbit plane immediately senses the direct effect of gravity as well 
as the more latent indirect effect (effect on trajectory) that conventional tracking data relies upon 
exclusively. 

The orbit determination aspects of a gravity mapping mission places different requirements on arc 
length as well. Precise range change measurements between two satellites do not provide sufficient 
information to determine independently the orbits of the two satellites. As will be demonstrated below, 
these measurements are extremely sensitive to some components of the orbits of the two satellites and 
almost completely insensitive to others (relative to how well they can be determined from other tracking 
data types). If the initial epoch state vectors of the two satellites are determined with other tracking 
data types, it should be possible to use the range change measurements by themselves to refine only 
certain components of the orbits. This has the potential to allow shorter arcs. Orbit refinement does 
not place the same constraints on arc length as orbit determination. 

Even if the use of short arcs in a gravity mission is possible, the question remains if their use is 
desirable. As mentioned above, the effect of unmodeled forces places an upper limit on the length of 
axes. In a mission which relies on high precision tracking, the tolerance for unmodeled forces is espe- 
cially low. Accelerometry can reduce the level of unmodeled forces, but it remains to be seen how well 
accelerometry will mix with very precise tracking data. For example, small thrusting events may be 
inadequately sampled by accelerometers or there may be small data gaps or noise. These situations are 
most straightforwardly handled by analyzing short arcs that avoid periods of anomalous accelerations. 
Furthermore, it is possible that for some applications, short arcs are desirable. For example, short arcs 
may facilitate the independent analysis of regional data in the form of gravity anomaly blocks or some 
other local gravity parameterization. It is beyond the scope of this investigation to answer many of the 
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questions surrounding the desirability of short arc analysis. We intend to demonstrate the possibility 
of short arc (under 15 minute in length) analysis. We also explore the tradeoff between using short 
arcs and arcs of more conventional length (1 day) assuming that the longer arcs are possible within 
the constraints presented by a real gravity mapping mission. 

3. Spherical Coordinates for Initial Epoch State Parameters 

Orbit solutions often solve for the six cartesian components of a satellite’s state vector at an initial 
epoch. After each iteration of an orbit solution, the updated initial cartesian state vector is used to 
compute the trajectory at later epochs, which is accomplished by numerical integration of the cartesian 
components of the state vector. Even though the initial epoch state vector is required in cartesian form, 
it is sometimes useful to solve for an alternative representation of the initial state vector for example, 
osculating Kepler elements. Often it is easier to apply useful constraint equations that are applicable 
to a particular data type when a non-cartesian representation is used. After the alternative form of 
the state vector is updated, it is transformed to cartesian coordinates and the numerical integration 
proceeds. It is straightforward to convert an orbit solution from a cartesian state vector solution to 
some other representation. All that is required is the six by six matrix of the partial derivatives of the 
six initial cartesian state vector parameters with respect to the six alternative parameters. 

In the case of the two satellites of a gravity SST mission, the twelve parameters describing the 
two initial cartesian state vectors can be converted to the cartesian state vector of the baseline (the 
difference vector) between the two satellites and the cartesian state vector of the baseline midpoint (the 
average of the two satellite state vectors). These two vectors can be further transformed. In the case of 
the midpoint state vector, it is useful to convert the three position components to spherical coordinates: 
the latitude of the midpoint, the longitude of the midpoint and the radius of the midpoint. For the 
baseline vector, it is useful to imagine a local cartesian coordinate system centered at the midpoint 
of the baseline. The XY plane of the system is perpendicular to the position vector of the system 
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midpoint. The X axis is the local East vector of the local coordinate system. It is further useful to 
describe both the position and velocity components of the baseline state vector in spherical coordinates 
which are based in this local coordinate system. In each case, the vector is converted to magnitude, 
pitch (angle the vector makes with XY plane) and yaw (angle that the projection onto the XY plane 
makes with the X axis). The twelve new epoch state parameters are: 

R m Distance of baseline midpoint from Earth center of mass 
Geocentric latitude of baseline midpoint 
A m Longitude of baseline midpoint 

X m Inertial X component of baseline midpoint velocity vector 

Y m Inertial Y component of baseline midpoint velocity vector 

Z m Inertial Z component of baseline midpoint velocity vector 

Rb p Baseline length 

Pbp Baseline position pitch 

Yb p Baseline position yaw 

R bv Baseline velocity vector magnitude 

P bv Baseline velocity vector pitch 

Ybv Baseline velocity vector yaw 


We have implemented the ability to solve for the above twelve parameters in our orbit determination 
and geodetic parameter estimation software, GEODYN (Pavlis et al., 2001). This required only a new 
subroutine which transforms back and forth between a pair of cartesian state vectors and the above 
twelve parameters and which also computes the twelve by twelve matrix of the partial derivatives of the 
twelve initial cartesian state parameters with respect to the twelve spherical coordinate parameters. 

In the next sections we demonstrate that the above twelve parameters are promising for use in the 
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analysis of the intersatellite range change measurements when they are decoupled from other tracking 
data types that might be available. One reasonably expects that the intersatellite range change 
measurement is more sensitive to the baseline parameters than to the midpoint parameters. One 
also expects that the intersatellite range change measurement is fairly insensitive to the baseline yaw 
parameters. If the range change measurements are to be analyzed independently from other tracking 
data types, then it is important to study the sensitivity of the range change measurement to each of 
the above twelve parameters. The sensitivity of the measurement to each of the twelve parameters 
should be compared with how well each parameter is likely to be determined independently from other 
tracking data types available in the mission. That subject is explored in Section 5. 

4. Data Simulation and Assumptions for Data Reduction 

All of the studies in this paper are based on using simulated one-way intersatellite range-rate 
data with a counting interval of 1 second. The data were generated using the orbits of two co-orbiting 
satellites with characteristics shown in Table 1. The key points are that the satellites are at an altitude 
very close to 500 kilometers, have a very low eccentricity, are very nearly polar and are separated by 
about 200 kilometers. The second satellite’s initial elements are an exact copy of the first satellite’s 
initial elements 30 seconds later. The satellite orbits were generated with the EGM96 gravity field 
using all terms up through degree 120. Drag and solar radiation were not modeled in our simulations; 
it is assumed that accelerometry will sufficiently account for these forces. Accelerometry should also 
compensate for small thrusting events, but it is possible that thrusting can cause occasional problems, 
so our simulations do include thrusting. Each satellite was given a thrust every 30 minutes (a AH 
of about 0.5 mms" 1 ) but in a staggered manner so that the satellite-satellite system received a AH 
every 15 minutes. (This, in fact, is roughly comparable to the thrusting frequency currently occurring 
on the CHAMP satellite [Schwintzer, 2000].) Although it would be hoped that an accelerometer 
would accurately model thrusts, we wish to determine if problematic thrusts can be removed from 
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the analysis process altogether (assuming they occur no more frequently than every 15 minutes). We 
therefore investigate whether gravity signal can be recovered while analyzing range- rate data in arcs 
of slightly less than 15 minutes (between thrusts). Baselines are refined every 15 minutes using arcs 
that begin 40 seconds after a thrust and 20 seconds before a thrust (14 minute arcs). 

9 

A data point was created every 5 seconds for 30 days. Two versions of the data were created, one 
with no noise and the other with noise of 1 //ms -1 . The noiseless version was useful in the verification 
of the modifications made to the GEODYN software for the new epoch state parameterization. The 
noiseless version was also used in preliminary studies to determine range-rate data sensitivity to epoch 
state parameters in orbit refinement solutions. The lack of noise makes residual analysis easier. All 
gravity recovery simulations used the data with noise. 

Most of the studies in this paper involve reducing the simulated data using an a priori gravity field 
which is different from EGM96. A crucial question is to determine whether the intersatellite range- 
rate data can be used by themselves (after orbits have been precomputed with a variety of tracking 
types and then refined in certain components with the intersatellite measurements) to recover the 
coefficients of EGM96, starting from a different gravity field. The a priori gravity field is a “clone” 
of EGM96 through degree 70 and is zero above that. At degree 70 and below, the coefficients of the 
a priori clone field differ from EGM96 by about the standard errors of EGM96. While considerable 
evidence suggests that the EGM96 errors are realistic (Lemoine et ah, 1998a), the EGM96 model 
itself does differ from some other recent models by amounts exceeding these errors. Our a priori clone 
gravity model should therefore be considered as fairly “close” to EGM96, so if we are able to recover 
a field that is closer to EGM96 than the clone, then we will have shown that our technique is useful 
for recovery of small gravity signal below degree 70. Basing the test on an a priori field that is already 
close to EGM96 is a stringent test of sensitivity. Above degree 70, our test will be less severe since 
we are starting from zero (which is farther away from EGM96 than the standard errors). Even so, we 
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will be able to determine if signal above degree 70 can be recovered with our method of reduction. 

The clone field was also used to generate the starting trajectory that is used to obtain the a priori 
elements to begin the orbit refinement process. The intersatellite range-rate data are used in the 
refinement step to adjust only three or four of the twelve initial state parameters (as discussed below). 
We accept the other eight or nine elements from the starting trajectory without alteration. So, in 
some sense, the starting trajectory contains the bulk of the the information about the orbits that will 
be used to extract gravity signal. 

In practice, the starting trajectories will have to be generated with whatever a priori gravity field 
is available. Our starting trajectory was determined in a manner somewhat approaching “reduced 
dynamic” methods (e.g., Bertiger et al., 1994; Rowlands et al., 1997). The solution had the following 
characteristics: 

1. The clone gravity field was used. 

2. AV's were solved for at the appropriate times. 

3. Empirical, periodic once per revolution accelerations (phase and amplitude parameters) were 
solved for in both the along- and cross-track components. Phase and amplitude were estimated 
every 30 minutes. 

4. The simulated intersatellite range-rate measurements were used as a tracking data type. 

5. To simulate the strong geometrical constraints that GPS tracking would provide, the “truth” 
ephemeris was used as a tracking data type. 

6. Solutions used thirty hours of data, from 0 hours of one day to 6 hours of the next. 

The starting ephemeris produced from the above solution differs from the “truth position” in an RMS 
sense by 3 centimeters (total position). That is probably optimistic for what will actually be attainable 
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for the position components of a reduced dynamic trajectory, but the velocity components are probably 
more important for this study (Jekeli, 2000). The starting trajectory fits the intersatellite range-rate 
data at 12 pms -1 , but that is a bit misleading. Although the starting trajectory was produced by 
reduced dynamic techniques, it will be used to generate a priori initial state elements for the gravity 
analysis. That analysis cannot make use of the empirical accelerations from the reduced dynamic 
trajectory (they contain gravity information). The a priori elements for the gravity analysis are 
gleaned from the starting trajectory and produce fits of over 200 pms -1 (RMS) in a 14 minute arc 
when they are used without the benefit of the empirical accelerations that were part of the starting 
trajectory. Our intersatellite range-rate gravity analysis starts with input elements which need to be 
refined along with the gravity. 

Note that empirical accelerations solved for while determining the reduced-dynamic starting trajec- 
tory do not alias into the gravity analysis, since the accelerations are used only to determine accurate 
intial conditions for the second (gravity estimation) step. The jump from 12 to 200 pins -1 mentioned 
above is a result of removing empirical accelerations in the second step. The 200 pms 1 contains 
signal from initial state error as well as from gravity errors that had been soaked up by the empirical 
accelerations. 

5. Baseline Refinement 

Before attempting any large simulations to demonstrate the ability of our technique to recover 
gravity information, we performed some smaller simulations with the goal of understanding baseline 
refinement from intersatellite measurements. In particular, we wish to determine which parameters 
need be refined and which parameters can be taken from the reduced dynamic trajectory without 
alteration. 

Our first set of tests used the truth force model (EGM96 through degree 120) and the noiseless 
intersatellite range-rate data. These tests were performed on short arcs (14 minutes) between the 
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AV thrusting events described in the previous section. The goal of these tests is to find the best 
minimum set of initial state parameters to estimate so that the a posteriori range rate residuals are 
well below the 1 pms~ l level (if all of the elements are set to truth values, then the noiseless data 
should fit perfectly when EGM96 is used). Although the data used in these tests were noiseless, they 
were weighted in the least squares solution as if they had a standard deviation of 1 micron per second. 
This is noted so that the formal standard deviations of the adjusted parameters can be interpreted. 

The first runs adjusted a single parameter. As noted in the previous section, elements taken from 
the reduced dynamic trajectories produce residuals with an RMS of over 200 pms~ l . When only a 
single parameter is adjusted, the only parameter that could reduce the RMS residual to under 100 
pms~ l was the velocity pitch parameter. In fact, the velocity pitch parameter adjustment produced 
an RMS residual of under 10 /ims -1 . 

The second set of runs adjusted two parameters with velocity pitch always being one of the pair. 
While examining the choice of a second parameter it became clear that over short arcs, two param- 
eters are very correlated: Rb v (baseline velocity magnitude) and Rm (the distance of the baseline 
midpoint from the center of mass of the earth). The choice of either of these parameters as the second 
parameter to accompany velocity pitch produces almost identical residual patterns. When these two 
parameters ( R bv and R m ) are allowed to adjust simultaneously along with velocity pitch, the inverted 
normal matrix shows a correlation between R n and Rt JV of very nearly one. When either of these two 
parameters accompanies velocity pitch, the solution produces an RMS residual of less than 1 pm s 1 . 
We chose to use velocity magnitude because the adjustments in this component were always less the 
100 /ims" 1 and usually less than 20 /ims" 1 . This seemed more reasonable than the adjustments in 
the R m parameter (sometimes more than 30 centimeters), since our reduced dynamic trajectories were 
better than 10 centimeters radially. In general the Rm parameter should almost always be determined 
from a reduced dynamic trajectory to better than 10 centimeters, so over short arcs this parameter is 
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not likely to need adjustment. 

Even though the adjustment of velocity pitch and velocity magnitude bring the RMS residuals to 
under a micron per second, there are noticeable trends in the residuals, sometimes near the micron per 
second level. Because of this we made one more set of runs to search for a third parameter. We find 
that baseline position pitch works best. With position pitch, velocity magnitude and velocity pitch 
adjusting, the RMS residuals were well under 0.1 pms -1 . Therefore, in our short-arc gravity recovery 
experiments (next section) we adjusted these three parameters. It may be possible to avoid adjusting 
position pitch, especially if arcs of 10 minutes or less are attempted. 

Our final set of runs deal with longer arcs. In the next section we make two gravity field determina- 
tions. One determination uses 30 days of 14 minute arcs (2878 arcs). The other uses 30 arcs, each one 
day in length. We want to find the proper set of adjusting parameters for one-day arcs. In extending 
to 2-hour arcs we found that our “14 minute parametrization” (position pitch, velocity magnitude and 
velocity pitch) held up quite well, producing fits of less than 0.2 pms -1 . We also found that for arc 
lengths of 2 hours, the R m parameter is still highly correlated with baseline velocity magnitude. At an 
arc length of 12 hours the “14 minute parameterization” produces RMS residuals of close to 1 pms -1 . 
Also, at this arc length, the R m parameter is less correlated (still 0.999) with the velocity magnitude 
parameter. At 12 hours and above, 4 parameters (including Rm) can be sensibly adjusted. In these 
adjustments the formal standard deviation of is less than a centimeter. When R m is adjusted as 
the fourth parameter in a 12 hour arc, the RMS residual is reduced to below a 0.1 /ims" 1 . In a 24 
hour arc, the correlation between R m and velocity magnitude is reduced to 0.995. Our “long arc 
gravity analysis described in the next section simultaneously adjusts four arc parameters along with 
gravity coefficients. 

As noted above, the noiseless version of the data was used in these baseline refinement studies 
and the force model was set to “truth” values. So if enough initial state parameters are allowed to 
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adjust, the data will be fit perfectly. The gravity analysis described in the next section uses the initial 
state parameterizations that have just been described, but the clone gravity field is used to compute 
the trajectories (refine the trajectories before normal equations are generated) and the “noise added” 
version of the data is used. In this mode, when the three short-arc parameters are allowed to adjust, 
the arcs fit the the data between 7 and 25 pms -1 . The day-long arcs, with 4 parameters adjusting, 
have fits between 20 and 60 jj, ms -1 . 

6. Gravity Solutions 

The gravity solutions presented in this section are, like most simulations, somewhat optimistic. 
The simulations include no unmodeled effects other than random measurement noise and some initial 
satellite state error (which is left over from those initial satellite state elements that we leave unaltered 
from the reduced dynamic trajectory). In the long-arc analysis we modeled (without adjustment) the 
“truth” values of the AV thrust events, which implies perfectly performing accelerometers. Futher- 
more, in addition to some initial satellite state refinement, only gravity parameters are estimated 
from the simulated data; no attempt is made to estimate, for example, tidal and gravity parameters 
simultaneously. In fact, tides and other high-frequency atmospheric and oceanic mass motions pose 
serious aliasing problems for an SST mission because of the difficulty in modeling and removing the 
associated gravity effects at required accuracies (Zlotnicki et al., 2000; Verhagen et ah, 2000). Such 
problems are here ignored. 

The two solutions presented below, one comprising short arcs and the other comprising long arcs, 
are intended to be taken qualitatively. The differences between the short-arc solution and the long-arc 
solution are of particular interest, since they reveal how much information is potentially lost when 
short arcs are used and what can be gained by extending arc length (assuming the level of unmodeled 
forces does not preclude the use of long arcs). But the fact that it is possible to obtain sensible gravity 
solutions from arcs shorter than 15 minutes is significant in itself. 
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We estimated two gravity fields using the 30 days of simulated data described in Section 4. Each 
gravity field was estimated to degree 120 without any constraints. The method used in estimating the 
fields differed in the choice of arc length. The “short arc” field used 2878 arcs 14 minutes in length 
while the “long arc” field consisted of 30 arcs, each a day in length. The short-arc field estimated 
2878 x 3 = 8634 arc (orbit) parameters simultaneously with the 14337 gravity coefficients while the 
long-arc field estimated 30 x 4 = 120 arc parameters. The short-arc field discarded 1 minute of data 
around each AF (12 points) every 15 minutes, so the short-arc field uses approximately 7% fewer 
observations. 

The estimated gravity fields should be compared with EGM96, which was used to simulate the 
data and is therefore the ‘true’ field, and the EGM96 clone, which was used as the a priori field from 
which gravity normal equations and initial ephemerides were produced. The figures in this section 
which pertain to coefficient values show differences from EGM96. If our estimates were perfect, the 
estimated coefficient differences from EGM96 would be zero. Of course the differences are not zero, 
but they are much smaller than the differences between EGM96 and the EGM96 clone (see Figure 2 
below). 

The formal errors of our estimated coefficients cannot be directly compared to EGM96 errors. 
The standard errors of EGM96 are the result of a complex calibration of weights of the many data 
types used in its solution (Lerch, 1991). The formal errors of our two estimated fields are simply 
the diagonals of an inverted normal matrix having a single data type which had been assigned a 
standard deviation of 1 micron per second. The formal errors of our estimated gravity coefficients 
should therefore be interpreted only in a relative sense. Internally, they should be reliable for seeing 
which portions of the estimated gravity fields are more strongly or weakly determined. Externally, 
they should provide a good basis to compare two gravity fields that were estimated in a largely similar 
fashion. 
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Figure 1 shows the RMS differences of the prior field and the two estimated fields with respect to 
the ‘true’ EGM96. It shows immediately (and reassuringly) that both estimated fields are considerable 
improvements over the a priori field. Neither field appears to be impacted by the truncation of the 
prior field at degree 70; both show relatively smooth differences with respect to EGM96 through all 
degrees. More interestingly, Figure 1 emphasizes that an accurate gravity field can be estimated from 
short arcs. The short-arc gravity field is significantly better than the a priori field at every degree 
from degree 5 through degree 100. It is not surprising that a field that is sewn together only from 
suborbital arcs is weaker at the very lowest degrees. The comparison of the performance (by degree) 
of the short-arc field with the long-arc field is not unfavorable to the short-arc field from degree 30 
upwards. However, at about degree 100, the short-arc field stops outperforming the clone gravity 
field. This does not happen until about degree 110 for the long-arc field. When judging the relative 
performance of the two estimated fields, it should be remembered that our simulations assume that 
the accelerometer is working perfectly (no unmodeled forces), which is much more beneficial to the 
outcome of the long- arc field than the short-arc field. 

More detailed comparisons of the three fields are shown in Figure 2, while the estimated formal 
errors are shown in Figure 3. Over a wide region of degrees and orders both estimated fields show 
remarkable improvements relative to the EGM96 clone. Many coefficients are improved by two orders 
of magnitude. This is generally consistent with figures quoted in the 1997 National Research Council 
report (see their Figure 2.6), but it is more definitive since the NRC calculations were based on 
an analytic theory of Jekeli and Rapp (1980) which assumes isotropic data (including data at all 
inclinations) and ignores possible required arc parameters. 

In general, Figures 2 and 3 agree well and show where the solutions are strong and where they 
are weakened when arcs are shortened. For a given degree, both the short-arc and long-arc fields 
determine lower order coefficients more accurately than higher order coefficients. That trend is much 
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more pronounced in the short-arc field, especially for sectorials, which for some low degrees are actually 
slightly inferior to the clone field. Again, this is not surprising— one cannot expect an extremely short 
arc (roughly 1 /6 of a revolution) in a high inclination orbit to be sensitive to long- wavelength sectorial 
terms. 

Most remarkably, the long-arc and short-arc fields are very comparable at low orders, especially 
so for zonal coefficients. The RMS discrepancy over all degrees (2 through 120) between EGM96 
and the short-arc field for zonal coefficients is 1.9 x 10~ 10 . For the long arc field that discrepancy 
is 1.8 x 10 -10 . Figure 4 shows that this striking similarity holds for all zonal terms, save the few 
between degrees 2 and about 8. For degrees 10 through 40 both long- arc and short-arc zonals are 
two orders of magnitude (or more) more accurate than the clone model. They are nearly one order of 
magnitude more accurate at degrees 40 through about 100. The improvement ceases at degree 112. 
Determining zonal gravity coefficients has historically been problematic in satellite geodesy. Clearly, 
an SST mission — even one which for one reason or another is restricted to using very short arcs of 
data — is likely to yield significant advances. 

7. Summary 

We have demonstrated a promising technique for the analysis of low-low intersatellite range-rate 
data from a gravity mapping mission. The technique largely (but not completely) decouples the task 
of orbit determination from the task of extracting gravity information. This has several advantages: 

1. In the first step of the procedure orbits can be determined using all available tracking types and 
reduced dynamic techniques. Extraction of gravity information in the second step benefits from 
the use of empirical accelerations used in the first step without aliasing problems. 

2. In the second step of the procedure, only intersatellite data are used. Gravity information can 
be extracted without complex solutions involving multiple data types. 
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3. Because orbits are only refined (not fully determined) in the second step, short arc analysis is 
facilitated. 

The technique transforms the standard 12 initial epoch state vector parameters of the two satellites 
into spherical coordinates describing the baseline between the two satellites. We have performed an 
analysis to show which of these 12 parameters need to be estimated simultaneously with gravity 
coefficients. 

We have used this technique to estimate gravity fields from simulated data. Our study neglects 
many effects, but unlike many studies of gravity missions in the literature, we consider the need to 
estimate orbit parameters simultaneously with gravity coefficients. We have investigated the possibility 
of estimating a gravity field entirely from short arcs (14 minutes) of intersatellite range rate data and 
found that this indeed should be possible. We have also shown the differences between a short-arc 
gravity field and a field estimated from long arcs. The use of long arcs (if possible) adds information 
primarily at the higher orders of every degree of the gravity field. 
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Table 1: Initial Keplerian Elements for Both Satellites 


Semimajor axis a 

6878050 m 

Eccentricity e 

0.001 

Inclination i 

89.1° 

Node n 

0° 

Perigee argument w 

0° 

Mean anomaly m 

0° 

Period T 

94.6 min 


Note: Initial epoch satellite 2 = initial epoch satellite 1 + 30 seconds. 


21 


FIGURE CAPTIONS 


Figure 1 . RMS gravity coefficient differences with respect to a “true” gravity field of (solid line) 
the a priori field, (dashed line) a solution employing short arcs, and (dotted line) a solution employing 
long arcs. The discontinuity in the a priori model arises because it is truncated to zero above degree 70. 


Figure 2. Results of gravity inversion simulations showing log y/(AC^ m + AS;( m ), for fully normal- 
ized Stokes coefficients C nm ,S nm differenced with the “true” gravity field, for (a) the a priori gravity 
field, (b) the short-arc gravity inversion, and (c) the long-arc gravity inversion. The discontinuity 
in the a priori model arises because it is truncated to zero above degree 70. Below degree 70 (a) is 
indicative of present-day uncertainties in the geopotential, as represented by the standard errors of 
EGM96. 


Figure 3. Formal errors for gravity simulations employing (a) short arcs and (b) long arcs. Colors 
show the logarithm of errors for fully normalized coefficients. 


Figure 4. Differences with respect to a “true” gravity field of zonal coefficients J n of degree n, 
for (solid line) the a priori field, (dashed line) the short-arc gravity solution, and (dotted line) the 
long-arc gravity solution. All three lines have been smoothed to remove minor statistical variability. 
The long-arc and short-arc solutions are similar except for the lowest degrees where the short-arc 
solution is less accurate. 
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